bbmleで最尤法

bbmleパッケージはRで最尤法を実行するためのパッケージであり、 optimなどの最適化ツールのwrapper関数であるmle2を使う。

データ例

ポアソン分布に従うカウントデータの例。 \[ y_i \sim {\rm Poisson}\left(\exp(\beta_0+\beta_1 \times x_i)\right), \quad i=1,\ldots,n \]

set.seed(100)
n <- 100
x <- runif(n,0,3)
y <- rpois(n=100,lambda=exp(0.1+1*x))
plot(x,y)
lines(sort(x),exp(0.1+1*sort(x)),lwd=2)

使用例

“\(-1 \times\) 対数尤度”に相等する関数を定義し、mle2関数で最適化（最小化）する。パラメーターの初期値startは名前付きリストで指定する。

minus.log.lik <- function(beta0,beta1){
  -sum(dpois(y,lambda=exp(beta0+beta1*x),log=TRUE))
}
library(bbmle)
res <- mle2(minus.log.lik,start=list(beta0=0,beta1=0))

summary関数，coef関数，vcov関数がすぐに使える。

summary(res)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = minus.log.lik, start = list(beta0 = 0, beta1 = 0))
## 
## Coefficients:
##       Estimate Std. Error z value  Pr(z)    
## beta0 0.157192   0.124194  1.2657 0.2056    
## beta1 0.935193   0.056309 16.6082 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 452.3947

coef(res)

##     beta0     beta1 
## 0.1571916 0.9351931

vcov(res)

##              beta0        beta1
## beta0  0.015424066 -0.006636444
## beta1 -0.006636444  0.003170724